R Notebook

Group 10

Libraries needed

Biological question and Project Scope

Clustering

kuang <- read.csv("./datasets/kuang-2014-microarray-expression.csv") ##local- needs to be changed 
time <- read.csv("./datasets/kuang-2014-microarry-expression-timepts.csv")
head(kuang[0:5])  ##preview data sets

tail(kuang[0:5]) 

rnaseq <- read.csv("./datasets/kuang-2014-rnaseq-expression.csv") #local needs to be changed
rnatime <- read.csv("./datasets/kuang-2014-rnaseq-expression-timepts.csv")

dim(kuang)

[1] 6241   37

# first remember the names
n <- kuang$Gene.name
# transpose all but the first column (Gene.name)
kuang_t <- as.data.frame(t(kuang[,-1]))
colnames(kuang_t) <- n
dim(kuang_t)

[1]   36 6241

dim(time)

[1] 36  1

kuang_t.time <- cbind(kuang_t, time = time$minute)
head(kuang_t.time[6241:6242])

colnames(kuang_t) %>% length()

[1] 6241

unique.genes <- colnames(kuang_t) %>% unique()
length(unique.genes)

[1] 6035

While looking at the data set we observed multiple genes that were present more than once. Using the unique function we can see that there are 6241 genes in the data set, but only 6035 unique genes. We are not sure for the reasoning behind this, our best guess is redundancy, but it is something to keep in mind.

Time plots

#visual of first 10 genes
first.10.genes <- c("YAL001C", "YAL002W", "YAL003W", "YAL004W", "YAL005C", "YAL007C", "YAL008W", "YAL009W", "YAL010C", "YAL011W")
ggplot(kuang_t.time, aes(x = time, y = YAL001C, group = 1)) + geom_line() + labs(x = "Time (minutes)", y = "Expression of Gene", title = "YAL001C")

ggplot(kuang_t.time, aes(x = time, y = YAL002W, group = 1)) + geom_line() + labs(x = "Time (minutes)", y = "Expression of Gene", title = "YAL002W")

ggplot(kuang_t.time, aes(x = time, y = YAL003W, group = 1)) + geom_line() + labs(x = "Time (minutes)", y = "Expression of Gene", title = "YAL003W")

ggplot(kuang_t.time, aes(x = time, y = YAL004W, group = 1)) + geom_line() + labs(x = "Time (minutes)", y = "Expression of Gene", title = "YAL004W")

ggplot(kuang_t.time, aes(x = time, y = YAL005C, group = 1)) + geom_line() + labs(x = "Time (minutes)", y = "Expression of Gene", title = "YAL005C")

ggplot(kuang_t.time, aes(x = time, y = YAL007C, group = 1)) + geom_line() + labs(x = "Time (minutes)", y = "Expression of Gene", title = "YAL007C")

ggplot(kuang_t.time, aes(x = time, y = YAL008W, group = 1)) + geom_line() + labs(x = "Time (minutes)", y = "Expression of Gene", title = "YAL008W")

ggplot(kuang_t.time, aes(x = time, y = YAL009W, group = 1)) + geom_line() + labs(x = "Time (minutes)", y = "Expression of Gene", title = "YAL009W")

ggplot(kuang_t.time, aes(x = time, y = YAL010C, group = 1)) + geom_line() + labs(x = "Time (minutes)", y = "Expression of Gene", title = "YAL010C")

ggplot(kuang_t.time, aes(x = time, y = YAL011W, group = 1)) + geom_line() + labs(x = "Time (minutes)", y = "Expression of Gene", title = "YAL011W")

Clustering

Complete method, with 1-Pearson correlation as the distance matrix

First, we attempted hierarchical clustering using the complete method with 1-Pearson’s correlation as the distance metric. The first thing we did was create the correlation matrix and then distance matrix of the transpose of the original data (with genes and the columns). We are clustered the transpose of the original matrix because we wanted to cluster genes based on their relative correlation over time. All of the clustering techniques we know cluster the columns, not rows, of the dataset.

## These take a lot of time so they are moved into their own code segment
kuang_t.cor <- kuang_t %>% cor(use="pairwise.complete.obs")  ##generate correlation measurements
kuang.dist <- as.dist(1 - kuang_t.cor)  ##use correlation to find distance measurements 

kuang.tree <- hclust(kuang.dist, method="complete")
kuang.dend <- as.dendrogram(kuang.tree) ## created & priviewed dendrogram object

To get an idea of the data we plotted the dendrogram. From the dendrogram it is clear that there are significant clusters in the dataset.

plot(kuang.dend, leaflab = "none")  ## removing labels 

To further get a feel for the dataset we plotted the dendrogram at several different cluster values (2, 4, and 8)- color coded to best analyze the data.

clusters_2 <- cutree(kuang.dend, k=2)
plot(color_branches(kuang.dend, k=2),leaflab="none")

table(clusters_2)

clusters_2
   1    2 
2956 3285 

clusters_4 <- cutree(kuang.dend, k=4)
plot(color_branches(kuang.dend, k=4),leaflab="none")

table(clusters_4)

clusters_4
   1    2    3    4 
2253  703 2313  972 

clusters_8 <- cutree(kuang.dend, k=8)
plot(color_branches(kuang.dend, k=8),leaflab="none")

table(clusters_8)

clusters_8
   1    2    3    4    5    6    7    8 
1535  703  404  718 1577  764  208  332 

sub.trees <- cut(kuang.dend, h = 1.75)  ##Height retrieved from 8 cluster Figure above
cluster1.final <- sub.trees$lower[[1]]
cluster2.final <- sub.trees$lower[[2]]
cluster3.final <- sub.trees$lower[[3]]
cluster4.final <- sub.trees$lower[[4]]
cluster5.final <- sub.trees$lower[[5]]
cluster6.final <- sub.trees$lower[[6]]
cluster7.final <- sub.trees$lower[[7]]
cluster8.final <- sub.trees$lower[[8]]
nleaves(cluster1.final)

[1] 703

nleaves(cluster2.final)

[1] 718

nleaves(cluster3.final)

[1] 1535

nleaves(cluster4.final)

[1] 764

nleaves(cluster5.final)

[1] 208

nleaves(cluster6.final)

[1] 332

nleaves(cluster7.final)

[1] 404

nleaves(cluster8.final)

[1] 1577

After reading our paper, we saw that they found three major clusters. We next did hierarchical clustering using k = 3 to see if our clusters align with the clusters found in the paper.

clusters_3 <- cutree(kuang.dend, k=3)
plot(color_branches(kuang.dend, k=3),leaflab="none", main = "Original Microarray dendrogram with k = 3")

table(clusters_3)

clusters_3
   1    2    3 
2956 2313  972 

We then cut the tree to extract these three clusters.

cut1 <- cut(kuang.dend, h = 1.9)
cut1$lower

Perfect. The first cluster countains 2956 genes, the second contains 972 genes, and the third contains 2313 genes.

To get a better idea we performed hierarchical clustering, starting with the first cluster on the left (pink)

plot(cut1$lower[[1]], leaflab="none", main = "Cluster 1 dendrogram")

We then tried plotting this at dendrogram at different k values to see what looks good.

plot(color_branches(cut1$lower[[1]], k=3), leaflab="none", main = "Cluster 1 dendrogram with k = 3")

The first two clusters (pink and green) looked pretty compact to us, so we did not explore them anymore. The third (blue) cluster appeared to have some significant subclusters in it, so we explored it more.

cluster1.cut1 <- cut(cut1$lower[[1]], h = 1.75)
cluster1.cut1$lower

[[1]]
'dendrogram' with 2 branches and 703 members total, at height 1.6383 

[[2]]
'dendrogram' with 2 branches and 718 members total, at height 1.708666 

[[3]]
'dendrogram' with 2 branches and 1535 members total, at height 1.7321 

plot(cluster1.cut1$lower[[3]], leaflab = "none")

We saw that there appeared to be 4 main clusters. Adding color allowed us to visualize if our observation held true.

plot(color_branches(cluster1.cut1$lower[[3]], k=4), leaflab="none", main = "Cluster 1 subtree dendrogram with k = 4")

That looked good to us. We then extracted these four clusters to see how large they were and what genes they contained.

cluster1.cut1.cutAgain <- cut(cluster1.cut1$lower[[3]], h = 1.4)
cluster1.cut1.cutAgain$lower

[[1]]
'dendrogram' with 2 branches and 312 members total, at height 1.361629 

[[2]]
'dendrogram' with 2 branches and 479 members total, at height 1.360901 

[[3]]
'dendrogram' with 2 branches and 342 members total, at height 1.138894 

[[4]]
'dendrogram' with 2 branches and 402 members total, at height 1.306937 

Next we looked at the smallest cluster ([[2]]).

plot(color_branches(cut1$lower[[2]], k=2), leaflab="none", main = "Cluster 2 dendrogram with k = 2")

Looking at this we saw there were two main clusters. We then cut this subtree again in order to extract these two subclusters.

cut1.cut2 <- cut(cut1$lower[[2]], h = 1.6)
cut1.cut2$lower

[[1]]
'dendrogram' with 2 branches and 764 members total, at height 1.413171 

[[2]]
'dendrogram' with 2 branches and 208 members total, at height 1.508391 

The smaller cluster (green) had 208 genes in it. This is a good size to work with, so we continued exploring this subcluster. First we plotted it.

plot(cut1.cut2$lower[[2]], leaflab = "none", main = "Subcluster of Cluster 2")

4 main clusters stuck out to us when looking at the above plot, so we plotted with color and k = 4.

plot(color_branches(cut1.cut2$lower[[2]], k=4), leaflab="none", main = "Subcluster of Cluster 2 dendrogram with k = 4")

We then extracted these four clusters.

cut1.cut2.cutTo4 <- cut(cut1.cut2$lower[[2]], h = 1.25)
cluster.genes <- as.character(labels(cut1.cut2.cutTo4$lower[[1]]))
cat(cluster.genes)

YOL146W YPL234C YDL097C YDR457W YDL224C YDR105C YFL020C YER152C YFR044C YKL135C YLL020C YLL020C YLL019C YLL020C YJL081C YJR017C YKL104C YAR019C YJL117W YBR016W YBR050C YDL159W YAR020C YCR038C YPL005W YPR091C YDR305C YDR366C YCR018C-A YLL061W YNL041C YNL140C YDR071C YPR188C YPR156C YCR097W-A YDR458C YOL109W YOL143C YOR122C YPL231W YLR447C YOL108C

kuang_t.cor.few <- kuang_t.cor[cluster.genes,cluster.genes]
f <- upper.tri(kuang_t.cor.few, diag=FALSE)
in.vector.form <- data.frame(vals = as.vector(kuang_t.cor.few[f]))
ggplot(in.vector.form, aes(x=vals)) + geom_histogram(bins=50)

The sizes looked reasonable to work with.

Lastly we analyzed the third ([[3]]) cluster (blue in the original dendrogram).

plot(cut1$lower[[3]], leaflab="none", main = "Cluster 3 dendrogram")

On this dendrogram we tried plotting it with different k values to see which made the most sense.

plot(color_branches(cut1$lower[[3]], k=5), leaflab="none", main = "Cluster 3 dendrogram with k = 5")

A k value of 5 looked good. As we have before, we then extracted these clusters by cutting the dendrogram.

cluster3.cut1 <- cut(cut1$lower[[3]], h = 1.6)
cluster3.cut1$lower

[[1]]
'dendrogram' with 2 branches and 332 members total, at height 1.565022 

[[2]]
'dendrogram' with 2 branches and 404 members total, at height 1.519815 

[[3]]
'dendrogram' with 2 branches and 141 members total, at height 1.563054 

[[4]]
'dendrogram' with 2 branches and 904 members total, at height 1.566672 

[[5]]
'dendrogram' with 2 branches and 532 members total, at height 1.5888 

RNAseq clustering

# first remember the names
nrna <- rnaseq$Gene.name
# transpose all but the first column (Gene.name)
rnaseq_t_zeros <- as.data.frame(t(rnaseq[,-1]))
colnames(rnaseq_t_zeros) <- nrna
#dim(rnaseq_t_zeros)
#dim(rnatime)
#tail(rnaseq_t_zeros[0:5])
rnaseq_t_zeros.time <- cbind(rnaseq_t_zeros, time = rnatime$hour)
#head(rnaseq_t_zeros.time[6241:6242])

Time plots for rnaseq

#visual of first 10 genes
first.10.genes <- c("YAL001C", "YAL002W", "YAL003W", "YAL004W", "YAL005C", "YAL007C", "YAL008W", "YAL009W", "YAL010C", "YAL011W")
ggplot(rnaseq_t_zeros.time, aes(x = time, y = YAL001C, group = 1)) + geom_line() + labs(x = "Time (minutes)", y = "Expression of Gene", title = "YAL001C")

ggplot(rnaseq_t_zeros.time, aes(x = time, y = YAL002W, group = 1)) + geom_line() + labs(x = "Time (minutes)", y = "Expression of Gene", title = "YAL002W")

ggplot(rnaseq_t_zeros.time, aes(x = time, y = YAL003W, group = 1)) + geom_line() + labs(x = "Time (minutes)", y = "Expression of Gene", title = "YAL003W")

ggplot(rnaseq_t_zeros.time, aes(x = time, y = YAL004W, group = 1)) + geom_line() + labs(x = "Time (minutes)", y = "Expression of Gene", title = "YAL004W")

ggplot(rnaseq_t_zeros.time, aes(x = time, y = YAL005C, group = 1)) + geom_line() + labs(x = "Time (minutes)", y = "Expression of Gene", title = "YAL005C")

ggplot(rnaseq_t_zeros.time, aes(x = time, y = YAL007C, group = 1)) + geom_line() + labs(x = "Time (minutes)", y = "Expression of Gene", title = "YAL007C")

ggplot(rnaseq_t_zeros.time, aes(x = time, y = YAL008W, group = 1)) + geom_line() + labs(x = "Time (minutes)", y = "Expression of Gene", title = "YAL008W")

ggplot(rnaseq_t_zeros.time, aes(x = time, y = YAL009W, group = 1)) + geom_line() + labs(x = "Time (minutes)", y = "Expression of Gene", title = "YAL009W")

ggplot(rnaseq_t_zeros.time, aes(x = time, y = YAL010C, group = 1)) + geom_line() + labs(x = "Time (minutes)", y = "Expression of Gene", title = "YAL010C")

ggplot(rnaseq_t_zeros.time, aes(x = time, y = YAL011W, group = 1)) + geom_line() + labs(x = "Time (minutes)", y = "Expression of Gene", title = "YAL011W")

dim(rnaseq_t_zeros)

[1]   16 6241

head(rnaseq_t_zeros[,6240:6241])

Some of these genes have zero expression and need to be removed!

head(rnaseq_t_zeros[, colSums(rnaseq_t_zeros == 0) > 0])

rnaseq_t <- rnaseq_t_zeros[, colSums(rnaseq_t_zeros != 0) > 0]
head(rnaseq_t)

rnaseq.cor <- cor(rnaseq_t, use = "pairwise.complete.obs") 
rnaseq.dist <- as.dist(1-rnaseq.cor)

rnaseq.tree <- hclust(rnaseq.dist, method = "complete") 
rnaseq.dend <- as.dendrogram(rnaseq.tree)

plot(rnaseq.dend, leaflab = "none", main = "Original RNAdeq dendrogram")

Compared to the microarray dendrogram, it looked like these clusters are more compact and therefore stronger. We predicted this based off of the lower heights of most clusters. In order to see how true this was, we explored the clusters more.

plot(color_branches(rnaseq.dend, k=3),leaflab="none", main = "Original RNAseq dendrogram with k = 3")